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Trajectory surface hopping (TSH) is one of the most widely used quantum-classical 
algorithms for nonadiabatic molecular dynamics. Despite its empirical effectiveness 
and popularity, a rigorous derivation of TSH as the classical limit of a combined 
quantum electron-nuclear dynamics is still missing. In this work we aim to elucidate 
the theoretical basis for the widely used hopping rules. Naturally, we concentrate 
thereby on the formal aspects of the TSH. 

Using a Gaussian wave packet limit, we derive the transition rates governing the 
hopping process at a simple avoided level crossing. In this derivation, which gives 
insight into the physics underlying the hopping process, some essential features of the 
standard TSH algorithm are retrieved, namely i) non-zero electronic transition rate 
("hopping probability") at avoided crossings; ii) rescaling of the nuclear velocities 
to conserve total energy; iii) electronic transition rates linear in the nonadiabatic 
coupling vectors. The well-known Landau-Zener model is then used for illustration. 
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I. INTRODUCTION 



Practical molecular modeling is based on the Born-Oppenheimer approximation (BOA), 
which allows one to decouple the fast electronic motion from the usually slower nuclear 
motion by introducing i) the (adiabatic) potential energy surfaces (PESs), provided by the 
electrons in a specific eigenstate, and ii) the nonadiabatic couplings (NACs).^ When the 
latter are neglected, nuclei move on a single PES. Despite the success of the BOA, there 
are many physical situations, such as photo-reaction, electron transfer, or any form of non- 
radiative electronic relaxation, which involve more than one PES.- In these cases one has to 
take into account the coupling between various PESs. 

One of the most widely applied techniques to treat nonadiabatic effects in molecular 
dynamics is the trajectory surface hopping (TSH), with its several variants.-"— The main idea 
behind this technique is that, while the electronic wave function is propagated coherently, 
the force field felt by the nuclei varies in a discontinuous, stochastic way — the nuclei 
move along a single adiabatic PES, selected according to the electronic population of the 
corresponding state; time changes in the electronic populations can result in a sudden hop 
to another adiabatic energy surface. In order to conserve the total energy after each hop, 
the nuclear velocities are rescaled. This leads to a discontinuity in the nuclear velocities 
which, however, is generally small since the hops are more likely to occur between PESs 
which are close in energy. Yet, energy conservation is obtained in a rather ad hoc way 
and, although it is common practice to rescale the velocities along the direction of the 
nonadiabatic coupling vectors (which couple different adiabatic states),-*^ in principle other 
choices are possible.— ^i^i^^ As a consequence, most surface hopping algorithms are justified 
on an empirical basis, by direct comparison with exact analytical results in model systems 
or experimental data. A further issue concerns the loss of electron-nuclear coherence in 
the course of the dynamics. This aspect is directly related to the scaling of the transition 
probability with respect to the nonadiabatic couplings. It is traditionally linear in standard 
TSH,"^ but it was recently shown to be incorrect for some model cases,— the effect being 
directly traced back to the treatment of decoherence over time. 

In this work emphasis is put on the formal analysis in order to provide a basis for better 
understanding surface hopping techniques, rather than concentrating on numerical aspects. 
This aim is similar in spirit to that of previous works,-'^^"— but in a simpler framework. 
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The paper is organized as follows. In Sec. Ullwe briefly introduce the formalism which will 
be used in the rest of the paper. In Sec. Illll we derive the equations governing the electronic 
transition rates at an avoided crossing by using a Gaussian wave packet limit for the nuclear 
wave function. The equations directly display the physics behind the TSH, i.e., the physics 
governing a hop at an avoided crossing followed by rescaling of the nuclear velocities in order 
to conserve total energy. We find that the physical source of the velocity rescaling is related 
to the speed of variation of the NACs and that the rescaling only affects the nuclear velocity 
components parallel to the NAG vectors. We further discuss some general consequences of 
our theoretical approach and its relevance in the context of practical nonadiabatic molecular 
dynamics. In Sec IIVI the electronic transition rates are explored by using the Landau-Zener 
model as a paradigmatic test case. We finally draw our conclusions and future perspectives. 



We consider a quantum mechanical system of n electrons and nuclei with the total 
Hamiltonian H = T + H being the sum of the kinetic energy of the nuclei, T = - Vr/2M, 
and the electronic Hamiltonian, H'^^, which contains the kinetic energy of the electrons, the 
electron-electron potential, the electron-nuclear coupling, and the nucleus-nucleus potential. 
To keep notations simple, throughout this paper we denote the set of electronic coordinates 
by r and the nuclear ones by R; moreover we consider all nuclei having the same mass M 
and we use atomic units. 

The nuclei are much heavier than the electrons and thus it is a natural first approximation 
to consider them as classical particles, having at all times well defined positions R and 
momenta P. This suggests an adiabatic procedure where the electronic problem is solved 
for nuclei momentarily clamped to fixed positions in space: H'^^(Il)(f>i{r; R) = i?i(R)(/)j(r; R). 
The (adiabatic) electronic eigenfunctions {0j(r;R)} depend parametrically on the atomic 
positions and form a complete and ortho normal set. They can be used as a basis to expand 
the total wave function of the system as 



II. BASICS 




(1) 



and solve the time-dependent Schrodinger equation 




(2) 
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The expansion coefficients Xji^] t): which depend on the nuclear positions, will be identified 
as nuclear wave packets. Such wave packets are neither orthogonal nor normalized. In fact, 
the integral f d^^R |xj(R'; = HXjli^)!!^ gives the instantaneous electronic population of 
the j*'^ quantum state. By inserting expansion ([T]) for the total wave function into the 
Schrodinger equation ([5]) and projecting out the resulting equation on the electronic state 
(/),•, one obtains 



-i5^D,,(R)-Px,(R;t). (3) 



with 

D^i(R) = ^(0.(r;R)|VR|0,(r;R)) (4) 

being the NAC vectors and P = — IVr the momentum operator for the nuclei, and where 
we used the notation (0j(r; R)|O|0j(r; R)) = J d^"r 0*(r; R)O0j(r; R), with O a generic 
operator. Note that in ([3]) we neglected the terms (0j(r; R)|T|0j(r; R)), being of the order 
1/M smaller than the kinetic energy of the electrons.— 

The NAC vectors Djj(R) couple different adiabatic energy surfaces. Generally the cou- 
pling terms in Eq. are of order 1 / a/M smaller than the electronic energy,— and therefore 
they can be safely neglected. In this case the adiabatic nuclear wave functions are evolved 
independently, i.e., their normalizations — which give the adiabatic populations — are 
constants of motion. The Born-Oppenheimer approximation is based on this decoupling. 
However, as the gap between two PESs narrows, the NAC vectors become large, as it can 
be seen from the following expression:— 

n 1 (0.(r;R)|(VRg°')|0,(r;R)) 

M E,(R)-E,(R) ' 

where -E'j(R) — -E'j(R) 7^ and the numerator remains finite. This allows mixing between 
eigenstates for large enough nuclear velocities. TSH then provides a convenient approxi- 
mation of the nonadiabatic molecular dynamics, i.e., of the electronic transitions that can 
occur along with the nuclear motion. 



III. BEYOND THE BORN-OPPENHEIMER APPROXIMATION 



A. Time-dependent perturbation theory 



A statistical reduction of a correlated electron-nuclear dynamics into occasional, indepen- 
dent, PES hopping is possible only if one transition (hop) is fully completed before the next 
can take place. This requires that the nonadiabatic coupling terms Djj(R), although not 
negligible, are small enough to be considered as a perturbation during a short time interval 



At, as typically assumed in analogous analyses, see, e.g., Refs. 



11, 



21 



and 
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In such a case, 

we can separate the Hamiltonian of the system in an unperturbed part H^^i = T + -E'i(R) 
and a perturbation part —^Ylj ■Dij(R) ■ P. The aim is to find an approximate solution of 
the time-dependent Schrodinger equation according to the time-dependent perturbation 
theory. In order to achieve this aim, it is convenient to work in the interaction picture^^, 
where 



X^,I = e''''-'x^ , 



Wiij = -i e 



Dy(R) P 



Note that in the adiabatic representation no time-ordering is needed in the definition of 
Wij. Therefore the solution of the equation of motion for Xi,i, between initial time zero and 
final time t, at first order in the perturbation reads 



x.,/(R;t) 0) 



- 5^ / dt' e'^"^'*'D,, (R) ■ P e-'^"-^*'xi(R; 0). 



(6) 



B. Approximations through Gaussians 

Eventually, ionic motion is treated classically while the computation of a hopping prob- 
ability has to proceed in a quantum mechanical framework. In order to establish the link 
between these two descriptions, we need a semi-classical approximation for the wave func- 
tions X- Inspired by Heller's work,— the initial state XjlR-jO) is represented in terms of 
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Gaussian wave packets 



^R,p,A,(R) = ( — ) exp 



A,- 



37V/4 



iP,-R-^(R-R,)' 



(7) 



as 

x,(R;0) = a,^R^,P,,A,(R), (8) 

with Rjo, Pjo, and aj being the average positions, average momenta, and amphtudes, re- 
spectively, of the wave packet at t = (see Appendix |X] for more details on the Gaussian 
wave packets). All Gaussians are normalized to one. The (complex) coefficient aj regu- 
lates the contribution from each Gaussian to the whole state. It expresses the correlations 
accumulated in previous time steps. The propagator e~^^°-^*' in ([6]) then evolves this wave 
packet from the initial time t = to a time t'. At this point we make use of a semi-classical 
approximation for the nuclei: in the spirit of the frozen Gaussian approximation proposed 
by Heller,— we impose that, for a short time interval, the width of the Gaussian is fixed 
("frozen") and that the time evolution of the parameters Rj(t) and Pj{t) is given by the 
solution of the classical equations of motion for an effective nuclear potential given by the 
jth pgg Note that the use of frozen Gaussians is a common practice both in numerical 
applications and formal developments in the field.— ^^i^i^"— One can then use the following 
approximation: 

e-^^°-*'6;R,oP,oA, ^ e-'(^^+^^)*'6;R^.(iOP,(t')A, , (9) 



p2 

with Tj = j^, Ej = Ejilij), and their sum Tj + Ej being constant along the classical 
evolution. In doing so, we are neglecting the term — dtPj(t) ■ Rj(t) in the quantum 
phase accumulated during the time evolution. This approximation is justified for small 
momentum changes during short-time propagation, which is in line with our derivation. In 
the following, we will refer to this semi-classical limit as the wave packet limit. Note that 
this limit is i n the sp irit of the short-time expansion to a semi-classical golden rule employed, 
271. and 



e.g., in Refs. 



25 



30 



For the sake of simplicity, we use a multivariate Gaussian with the same (frozen) width 
Xj for all nuclear Cartesian coordinates. This is justified in the classical limit Xj — )■ oo, 
which will be taken at the end of our derivation. Moreover in the following we will drop the 
time dependence of the average positions and momenta, if not needed. At each time the 



Gaussian wave packets fulfill the completeness relation 



where X^^\. is the inverse of the overlap Xn.,x,(Pk,Pe) = (^R^PfeA, I^r>p«a>) (see Appendix |X] 
for details). Inserting (jS]) and ( fTOl) in ([6]) and applying the wave packet limit ([9]) yields 

Xi,,(R;i)~x.(R;0)-5^a,^ dt' j d^^p^ j d=^^p, 

e'^^-^^^-^-^^-^*'^R.P.oA.(R)X^J,^(P.,P.) 

(^R,P,AjDiilfe,PjAj) • Pj , (11) 

where we applied the wave packet limit also to the P operator in ([6]) allowing the identifi- 
cation P = Pj. Note that, in Eq. (HH), ^r^oP^^a, = e-^(^»+^'=)*'e'^«--*'^R,(t/)p,(t.)A, implicitly 
depends on t' as initial time of the backward evolution of (Ri(t), Pfc(t)) from t' to t = 0. 

We have now to decide how to deal with the NAG vectors Dij. The adiabatic basis is 
associated with strongly varying Djj(R). Therefore, we will consider the following Gaussian 
distribution for the coupling vectors 

D,,- (R) = dJ^'^ exp [- (R - Rf )) ^ /i^'^) (R - R^'^^) 1 , (12) 

where T denotes transposition, Dq"*^ is a constant, and Rc*"''' the position of the avoided 
crossing. For notational convenience we will drop the superscript "(zj)" on the right-hand 
side of Eq. f lT2|) . The Gaussian "width" fi^'^^^ is a rank 2 tensor, whose form will be discussed 
in Sec. IIII G[ Ansatz (fT2|) is in line with the avoided crossing model proposed, e.g., in Ref. 



and with the analysis of Ref. |21| based on a semi-classical propagator; moreover it allows the 
NAGs to fulfill the curl condition,— at least in the case of a two- level system (2LS), as it is 
illustrated in Sec. IIII G[ 

The Gaussian form for D^- allows an analytical evaluation of the transition matrix el- 
ements. Moreover, to keep contact with the TSH technique, we consider that transitions 
j — )■ 2 at an avoided crossing produce again wave packets of about the same spatial width 
(i.e., Xi = Xj = X) and same average position (i.e., Rj = Hi). Therefore, using the folding 



relations of Gaussians and the inverse Xj^^^_)^(Pfc, P^) given in Appendix El one obtains 



Do ■ Pj exp 



-(P,-P,)V'(P, 



(13) 



^ a/ det (/i) Jo 

1 

In principle, one cannot move the exponentials out of the time integral because the wave 
packet parameters — being evolved according to the classical equations of motion — can 
display a non-trivial time dependence. On the other hand, we eventually consider the 
classical limit A —t- 00 of the previous equation. In this limit, one can consider the evolution 
of the wave packet parameters to be smooth over a time scale, t, large enough to approximate 
the time-integral with a Dirac delta-function. This is also justified in the proper classical 
limit because energy fluctuations are suppressed, as in classical molecular dynamics the total 
energy is exactly conserved at each time-step. The result then becomes 

x.(R;t)^x.(R;0)-(-) E^/^^P^ 

^R,oP.oA(R)e^^^^-^^^^^Do ■ P, 

^Ie^ + ^-E. '-]. (14) 



exp 



i(p,-p,)Tr'(p,-Pfc) 



2M ' 1M ^ 

We can now take the limit A — ?■ 00, which simply localizes the Gaussian wave packet 
^RioPfeoA(R) while leaving unchanged the rest of the expression. 

C. The curl condition for the NACs 

The NAG vectors are known to satisfy the so-called curl condition if they are not in the 
neighborhood of a conical intersection.^^ For an arbitrary number of electronic PESs, the 
curl condition is nonlinear in the components of the NAG vectors, and its analysis is beyond 
the scope of this article. However, for a 2LS with real wave functions the curl condition is 
linear, and reads 

aDP)(R) dDf\^) _^^ ^^^^ 

where Da'^\'R) are the components of the single nonzero independent NAG vector of the 
system, D*^^^)(R) = — D*^^^^(R), and the equation holds for all pairs of nuclear coordinates 
(a,/3). 
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Ansatz (fT2!) is flexible enough to adapt to the curl condition for a 2LS, Eq. (fT5|) . First, fi 
should be a semi-positive symmetric matrix in the nuclear coordinates, i.e., it should satisfy 
= l^pa and R^'^/tR > for all R. The properties of such a matrix fi guarantee that 
i) there is a change of nuclear coordinates associated to an orthonormal basis change (i.e., 
a rotation) which transforms fi into diagonal form, and that ii) some of its eigenvalues are 
positive, and some may be zero. Directions with zero eigenvalues of fi give rise to Dirac deltas 
in momentum space instead of flnite-width Gaussians, and hence there is strict momentum 
conservation along these directions; eigen-directions with nonzero eigenvalues are described 
via Gaussians in momentum space, with the fi restricted to this invertible subspace, and 
hence small changes of the nuclear momentum along these directions are allowed. 

It is possible to prove (see Appendix |B]) that for such an ansatz the curl condition is 
satisfled for all nuclear conflgurations R if and only if 



with a positive proportionality constant. Such a /i has a single nonzero eigenvalue, which 
corresponds to the direction of Dq. In the following we will consider a /i as given in ( IT6l) 
also for the general case of multiple PESs. 

D. Transition rates 

From the flnal result f ll4p of perturbation theory in the Gaussian wave packet approx- 
imation, we can derive the change in time of the electronic population in the i*^ state, 
llx«(^)lP ~ IIXi(0)lP? from which the electronic transition rates between an initial state Xj 
and the flnal states Xi are obtained as 



Equation ( |T71) is the central result of this paper: it describes hopping between an initial 
adiabatic energy surface Ej, along which nuclei move with momenta Pj, and a flnal adiabatic 
surface Ei, along which nuclei move with rescaled momenta in order to conserve the total 
energy. This is precisely the framework common to the various TSH approaches. 



/t oc Do ® Do 



(16) 



jP.^iPfc Re«aj) Do ■ Pj exp -^(Pj - ^kff^ ^(Pj - Pfe) 




(17) 
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Besides recovering the essential features of the TSH algorithm, our derivation provides a 
better understanding of the underlying physics. In particular, the change Pj — )■ in the 
nuclear momenta occurs within a range set by /t, which is related to the spatial variation 



21 



We note 



of the nonadiabatic coupling vector. A similar result can be found in Ref. 
that when considering nearly constant Dj^, i.e., when fi — > 0, the exponential in Eq. f ll7p 
becomes 6^^{Pk — Pj) and the energy matching becomes 6 {Ei — Ej) requiring a strict level 
crossing. This is the case in a diabatic basis, as we will see in Sec. IIVI 

The allowed changes in momentum are aligned along the direction of Dq, i.e., the direction 
of the NAC vectors, as discussed in Sec. IIII CI This result supports hence the widely used 
procedure of adjusting the nuclear velocities along the NAC vectors and it is in line with 
previous findings in this direction.— Note that this result strongly relies on the ansatz 
(fT2|) and (fT6|) for the NACs. This choice allows the NACs to satisfy the curl condition for 
a 2LS (with real wave functions), Eq. f fTSj) . and in our derivation we assume the same form 
also for a general multi-level system. 

Finally, we find a transition rate linear in the coupling vector Dq, typical of the standard 
surface hopping algorithm.^ Crucial in getting this scaling is to assume that the state Xi is 
initially populated in Eq. ^ which means a non-zero ai in Eq. (1171) . This requires that the 
electronic correlations contained in the aj coefficients were propagated coherently over the 
history of the process. 

If one requires, instead, that a full state reduction is performed at each time when evalu- 
ating the transition rates, then each coherent propagation starts from a pure single state j', 
i.e., aj = 6jj>, which removes the sum over j in f[T^ . In this case the transition rates to 
previously unoccupied levels i ^ j' are quadratic in Dq, since the linear term drops from 
Eq. ([IT]): 



W^j'p,.,^iPfc oc |Do ■ Pj'P exp 



-l{Py-P,ff,-\Py 



2 

p2 p2 

SlE, + ^-Ei, '-]. (18) 

^ ' 2M ^ 2M ' ^ ^ 



These results are in line with the recent findings that Tully's surface hopping gives the wrong 
scaling in the spin-boson model (the correct scaling being quadratic in the coupling vector) 
and that this is due to an incorrect description of decoherence in the standard TSH.— 
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IV. ILLUSTRATION IN THE LANDAU-ZENER MODEL 



The final formula ( [T7|) can be illustrated in the well known Landau-Zener model.— The 
basics of the model are sketched schematically in Fig. [TJ The nuclear degree of freedom is 
described by one coordinate R. The electronic degrees of freedom are focused to a system 
of two levels i = 1,2. The unperturbed system, standing for the diabatic situation, has a 
linear level crossing at i? = 0. The slope E' is the first model parameter. The strength 
of the interaction between the two diabatic levels is set by the coupling constant V, which 
constitutes the second model parameter. The adiabatic representation is obtained by solving 
the 2x2 model Hamiltonian for fixed ionic position. This yields the well known adiabatic 
PESs El 2 = ±a/ {RE'Y + V"^ as indicated in the figure (solid lines). These are the PESs 
which enter the hopping formula f|T7|) . The deviation from the diabatic energy levels is 
particularly strong around i? ~ 0, which leads to a strongly varying nonadiabatic coupling 
Djj(R) as was assumed in our derivation. It is a textbook exercise to work that out for the 
Landau-Zener model. We find 

This matrix element is strongly peaked at the avoided level crossing, i.e., around i? = 0, 
with a characteristic R width This translates to a typical width of the 

momentum distribution of The example demonstrates that some non-negligible 

momentum spread in the hopping is expected as we usually encounter avoided crossings 
as modeled in the Landau-Zener model. The overall strength of hopping matrix element 
is governed by the same parameter combination l-E'/^l which determines the momentum 
width. We thus find that larger hopping probabilities are associated with larger momentum 
widths. 

One can also try to generalize the Landau-Zener model to higher dimensions. A simplistic 
way to achieve this is to promote E' to a vector, and, consequently, to interpret RE' as a 
scalar product. In this case Eq. (1191) indicates that the direction of D12 is along E' and 
/X is proportional to the tensor product E' ® E' . This further supports ansatz (fT2l) . with 
fi cx Do ® Do, for the NACs, and our findings of Sec. IIII D[ 
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FIG. 1. Schematic plot of the diabatic energies zizRE' (dotted hnes) and adiabatic energies Ei^2 
(sohd hnes) in the Landau-Zener model as functions of the ionic distance R. 

V. CONCLUSIONS 

We propose a derivation of the trajectory surface-hopping (TSH) technique based on a 
semi-classical approximation to the nuclear dynamics in the spirit of a wave packet limit. 

The equations governing the electronic transition rates at a simple avoided level crossing 
display the essential features of TSH algorithm and allows us to elucidate the underlying 
physics. We find a nonzero electronic transition rate at avoided crossings, which allows for 
small changes in the nuclear momenta accounted for properly in the energy matching. This 
justifies the rescaling of the classical velocities done in practice after each hop. Moreover, 
we find that the physical source of the width of allowed hops in momentum space is related 
to the speed of variation of the nonadiabatic coupling elements in the adiabatic basis. In 
the classical limit for the nuclei the derivation supports the rescaling of the momenta along 
the nonadiabatic coupling vectors. This result strongly relies on the ansatz employed for 
the nonadiabatic couplings (NACs). In our derivation we assume a multivariate Gaussian 
form for the NACs, which allows the NACs to fulfill of the so-called curl condition, at least 
in a two-level case. We also find that the final electronic transition rate is linear in the 
nonadiabatic coupling vectors, as in the standard TSH algorithm, and that incorporating 
quantum decoherence makes this scaling quadratic. 

Illustration through the Landau-Zener model supports our findings. 
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Appendix A: Gaussian wave packets 

In this appendix, we collect a few general properties of the Gaussian wave packets, for 
which most integrals are known analytically. For example, the folding of Gaussians in general 
obeys the simple rule 



abnV^^^ fhR' + aK" 
exp iP 



a + b J \ a + b 

f oh . (R'-R")'\ , 
V 4(a + 6) a + b J ^ ^ 

The basic multivariate isotropic Gaussian wave function reads 

e^R.p,A(R)= - exp iP,-R--(R-R,)M . (A2) 



where A controls the spatial width of the Gaussian. In the present work we only consider 
overlaps between two Gaussian wave packets with the same spatial widths and centers, 

^ gi(p,-p.) R. ^_ . (A3) 

The inverse of these overlaps, Z^^^(Pf:,Pj) defined to satisfy 

I d3^PfcXa,,(P„ P.)X^;,(P,, P,) = 5=^^(P, - P,) , (A4) 
may be expressed as 

3N/2 

e 



/ 1 \' 

^K:.(p<.p.)=(sAj 



One can verify that such an inverse also satisfies the completeness relation (fTOlK 
Details on more general multivariate Gaussians can be found, e.g., in Ref. 
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Appendix B: Proof of the curl-consistency of jl for a two- level system 



In the basis where /i is diagonal, 

■iN 



D(R) = Doexp<^-5^/i4(R-Re)jn , (Bl) 



SO 

dDJ-R) dD.CR) 



-2 



/i^(R-R,LZ},(R)-/i,(R-Re),/^^(R) . (B2) 



If we enforce the curl condition for two levels at the point R = Rc + ce^, where (R — Rc)a = 
we get 

dp^fi^D^]^^ exp (-/i^c^) = 6a-yfi-yD^Q^f exp (-/i^c^) (B3) 
for all pairs of directions (a, f3). If we take /3 = 7, we have 

for all components a. There are only two ways to satisfy these equations: either /x^ = 0, 
or -Do!a^ = ^aj^o^,^^ ^ot- This implies that, if 7 is a direction such that /i^ 7^ 0, then, for all 
components a 7^ 7, ^ = 0. Since, in order to have nonzero Gaussian NAG vectors, at 
least one component of D and one eigenvalue of fi should be different from zero, there can 
only be one non-zero eigenvalue of fx, and it will correspond to the same direction of the single 
non-zero component of Dq . Therefore, in a base-independent expression, fi oc Dq ® Dq, 
where the tensor product is defined in terms of components as (a ®h)aii = a-abfi, and the 
proportionality constant must be a positive real number. It is straightforward that this 
condition is not only necessary but sufficient, since such a fi will always satisfy Eq. ( 1B2I) . 
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